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We carry out a layer-by-layer investigation to understand electron transport across metal- 
insulator-metal junctions. Interfacial structures of junctions were studied and characterized using 
first-principles density functional theory within the generalized gradient approximation. We found 
that as a function of the number of crystal layers the calculated transmission coefficients of multilayer 
silicene junctions decay much slower than for /i-BN-based junctions We revisited the semiclassical 
Boltzmann theory of electronic transport and applied to multilayer silicene and /i-BN-based junc¬ 
tions. The calculated resistance in the high-transmission regime is smaller than that provided by 
the Landauer formula. As the thickness of the barrier increases, results from the Boltzmann and the 
Landauer formulae converge. We provide a upper limit in the transmission coefficient below which, 
the Landauer method becomes valid. Quantitatively, when the transmission coefficient is lower than 
~ 0.05 per channel, the error introduced by the Landauer formula for calculating the resistance is 
negligible. In addition, we found that the resistance of a junction is not entirely determined by the 
averaged transmission, but also by the distribution of the transmission over the first Brillouin zone. 
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I. INTRODUCTION 

Silicene is a single atomic layer of Si atoms arranged in a two-dimensional buckled honeycomb lattice. Although it 
was predicted decades ago that free-standing silicene exhibits a massless relativistic behavior near the Fermi energy,— “— 
it was only very recently that silicene was synthesized on the surfaces of a few metallic substrates— Linear electronic 
dispersions observed in the silicene/Ag(l 11) system were attributed as the signature of Dirac fermions in silicene— 
Subsequent experimental and theoretical studies revealed the absence of Dirac fermions near the Fermi energy; instead, 
the observed linear electronic dispersions are sp-bands of Ag or hybrid interface states.— “— 

Currently there is a growing interest in multilayer silicenes.— The atomic structure of multilayer-silicene has 
not yet been identified experimentally. The surface morphology of multilayer-silicene on Ag(lll) was examined using 
scanning tunneling microscopy (STM) and scanning electron microscopy (SEM). Growth of the first silicene layer 
on an Ag(lll) substrate forms a 4 x 4 super cell with respect to the Ag(lll) surface and 3x3 with respect to 
free-standing silicene. For multilayer silicene, a (v^ x •\/3)R-30° reconstruction with respect to the free-standing 
silicene was observed— The reconstruction in multilayer-silicene/Ag(lll) was reproduced theoretically, and it was 
predicted that only the surface silicene layer reconstructs 

The electronic transport properties of silicenes are important for their possible applications. The current-in-plane 
(CIP) configuration is a natural choice to measure the resistance of silicenes, but the current tends to pass through the 
highly conducting Ag(lll) substrates, which makes measurements of the CIP resistances rather difficult—^ Moreover, 
multilayer silicenes grow in the so-called Stranski-Krastanov mode (also known as the layer-plus-island mode), in 
whicbi^d^ the first silicene layer forms a continuous film on Ag(lll), while subsequent Si atoms form islands of 
multilayer-silicene. As a result, it is difficult to directly associate the measured CIP resistance with the thickness of 
a multilayer silicene. It is much easier to perform a current-perpendicular-to-plane (CPP) measurement by placing a 
conducting probe on top of a multilayer-silicene island and measuring the resistance between the probe and the Ag 
substrate. 

The transport properties of ultrathin silicene-based junctions were studied previously,— and the calculated average 
transmission per channel are about 0.7 and 0.35 for monolayer and bilayer silicene based junctions, respectively. The 
Landauer formula is no longer accurate for calculating the resistances of junctions with such high transmissions, 
and we instead applied the semiclassical Boltzmann (SCB) theory to calculate the resistance. In this work we 
extended our previous work to study silicene-based junctions with up to eight silicene layers. The comparison between 
calculated SCB and Landauer resistances reveals the relation between them, and more importantly provides an 
empirical threshold for the application of the Landauer formula. 

The rest of the paper is organized as follows. The computation method is presented in Section [III The calculation 
results for Ag(III)|multilayer-silicene|Ag junctions are shown in Section lllll Complementary calculations on multilayer 
/i-BN-based junctions are presented in Sec. IIVI A summary is given in Section IVl 


II. COMPUTATIONAL METHOD 


The atomic model for Ag(III)|multilayer-silicene|Ag(lll) junctions was built according to the first-principles sim¬ 
ulations of bilayer-silicene/Ag(lll) interfaces presented in Ref. without considering atomic reconstructions at the 
interfaces. The junctions with two to eight silicene layers were optimized using the projector-augmented wav o^^i^^ 
(PAW) based density functional theory (DFT) as implemented in the Vienna ab initio simulation package VASP— 

In this work we employed the generalized gradient approximation (GCA) with the Perdew-Burke-Ernzerhof (PBE) 
parametrizatiorii^ The supercells used in the structural optimizations consist of the multilayer silicene and five Ag(III) 
atomic layers with Ag atoms in the central atomic layer kept fixed in their bulk positions. The lattice constant of 
free-standing monolayer-silicene is approximately 4/3 times of the Ag(III) surface, and a supercell consisting of 
4x4 Ag(lll) layers and 3x3 multilayer-silicene primitive unit cells in the x-y plane was used to simulate these 
junctions— The lattice constant in the x-y plane was fixed by that of Ag(III), calculated using the PBE functional 
to be 2.952 A; and the lattice constant of multilayer-silicene was 2.952 A x 4/3 = 3.936 A. During structure optimiza¬ 
tions, the height of the supercell (along the 2 -direction) and the coordinates of atoms were fully relaxed, until the 
forces on unfixed atoms were smaller than 0.01 eV/A. 

The DFT-based NEGFi^SA^S method was then used to compute the Green’s function and the total transmission of 
these junctions. NEGF calculations were performed using the transiesta code^. Numerical atomic orbitals were 
used to expand the Hamiltonian and the Green’s function. Single-zeta plus polarization orbital (SZP) and double-zeta 
plus polarization orbital basis sets (DZP) for Ag and Si respectively were generated using the default parameters in 
SIESTA— Norm-conserving pseudopotentials^ were used to describe interactions between valence electrons (3s^3p^ for 
Si and 4d^°5s^ for Ag) and the corresponding core electrons. The direction of transport was chosen as the z-direction. 
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FIG. 1. (Color online) Schematics of a system in the current-perpendicular-to-plane (CPP) configuration. The system consists 
of n metallic leads which are separated by n — 1 junctions. The system is homogeneous in the x-y plane, and the electrical 
current flows in the 2 direction. There are only two leads and one junction, i.e., n = 2 in our calculations on multilayer silicene 
and /i-BN-based junctions. 


Translational symmetry in the x-y plane was exploited by using a 15 x 15 fc-point mesh for calculating the charge 
density and 55 x 55 for the Green’s function. 

The group velocity of the Bloch states in the Ag(lll) leads along the 2 -direction is2^ 

ui(k||) = ^[u^(k|,)]^r'-V(k||), (1) 


where is the length of the unit cell of the Ag lead; U'^ (k||) is the periodic part of the Bloch waves; ky are the 
components of fc-points in the x-y plane; j is the index for Bloch waves at the same ky; and — E“), where 

S’' and are the retarded and advanced self-energy of the Ag leads, respectively. The transmission t and reflection 
r coefficients corresponding to the Bloch waves in the left and the right leads are^^ 
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where the dependence on ky is omitted. The label M denotes one of the Ag leads (left or right) and N the other one; 
the subscripts > and < denote Bloch waves propagating against and towards the junction, respectively; and is 

the submatrix of the retarded Green function of the junctions. 

The group velocities of Bloch waves in the leads [Eq. ([T|)] and the transmission and reflection coefficients [Eqs. ([2]), 
(|3|)] extracted from first-principles calculations are used as parameters for the Boltzmann equation. It is convenient 
to introduce an auxiliary quantity h to characterize the change of the distribution function / from its equilibrium /o 
with energy, 


/^(2,k||,£;)=/o(E)-/l^(2,k||)||. 


The Boltzmann equation for the CPP configuratioi 
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where r is the relaxation time in Ag leads; y,{z) and £z denote the applied electric field and the chemical potential 
along the 2 -direction, respectively. The method to solve the Boltzmann equation Eq. ([5]) numerically is given in the 
Appendix. 

The current density along the 2 -direction is 

'^ = -^^^sgn[ui(k||)] /i^'(2,k||) = ^Jk„(2), (6) 
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with A the cross section of the unit cell of junction perpendicular to the z-direction. The total current density J is a 
constant due to the conservation of charge, although each of its components Jk|| (z) is not necessarily a constant. The 
expression for the local chemical potential (see Appendix B) is ^(z) = {h^{z, k||))^. , and the voltage drop across the 

junction located at z = zq is equal to AV = ^(zq + 0+) — ^(zq — 0“); thus the four-probe resistance of the junction 
is calculated as AV/ J. 


III. MULTILAYER SILICENE JUNCTIONS 





AA-Si 


ABC-Si 


FIG. 2. (Color online) Atomic structures of bulk AA-Si and ABC-Si. The unit cell boundary is denoted by blue lines. 


We considered two different stacking orders for the multilayer-silicenesfi^ In the first stacking order, two inequivalent 
silicene layers are stacked in an “AA” manner (denoted as “AA-Si”), and each Si atom can find another in its 
neighbouring layers with the same in-plane position. The second stacking order, denoted as “ABC-Si”, corresponds to 
the stacking along the (111) direction of diamond-structured silicon. We note that both of these stacking configurations 
lead to a tetragonal arrangement of Si atoms, as shown in Fig. [21 
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FIG. 3. (Color online) Energetics of (a) bulk AA-Si and ABC-Si, and (b) Ag(lll)|AA-Si|Ag and Ag(lll)|ABC-Si|Ag junctions 
as a function of the number of Si layers. 

The total energies of bulk AA- and ABC-Si were calculated as a function of the interlayer distance. The interlayer 
distance is one half (third) of the lattice constant of bulk AA-Si (ABC-Si) along the z-direction. Bulk AA-Si has a 
higher total energy than bulk ABC-Si by 14meV per Si atom, as shown in Fig.|3Ka). 

We also calculated the total energies of multilayer silicene based junctions. The total energies of AA-Si based 
junctions are always higher than the corresponding ABC-Si junctions, and the total energy difference per Si atom 
is shown in Fig. |3l(b). The total energy difference between bulk AA-Si and ABC-Si is denoted as the dashed line in 
Fig- mb). The deviations from the dashed line in Fig.|3l(b) are due to the interface effect. We had also tried another 
method to obtain the atomic structures of junctions, in which the junction with {N + 1) silicene layers was constructed 
and optimized from the junction with N silicene layers by inserting a flat Si layer between the Nth layer of silicene 
and the Ag(lll) lead. The resulting multilayer silicene structures are different from either AA-Si or ABC-Si, but are 
similar to the body-centered tetragonal C 4 allotrope of carboii2£. Results of these structures are not shown due to 
their much higher energies. 
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A. Transmission 
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FIG. 4. (Color online) The averaged transmission over kx for Ag|AA-Si|Ag and Ag|ABC-Si|Ag junctions as a function of the 
number of Si layers. 

The transmission of Ag(lll)|multilayer-silicene|Ag(lll) junctions at the Fermi energy were calculated using the 
DFT-NEGF method, averaged over k||-points in the first Brillouin zone, 

r=^Er(k,). (7) 

" kll 

Because there is more than one transverse mode for each kx due to the band structure folding, the averaged trans¬ 
mission T could be larger than unity. For thinner junctions, with less than four silicene layers, the transmissions of 
AA-Si junctions are very close to those of ABC-Si junctions. For thicker barriers, the averaged transmission decays 
exponentially as a function of the number of silicene layers. The decay rate of transmission of ABC-Si junctions is 
smaller than that of AA-Si junctions. 



FIG. 5. (Color online) The density of k, denoted as n(K), for (a) AA-Si and (b) ABC-Si respectively; k is defined as 
K(k||) = mink|| \m{kz)\ and the relative importance of kI{k,) = n(K)e“”‘*/[n(Kinin)e“'‘“‘”‘*] for (c) AA-Si and (d) ABC-Si 
respectively, with d — 2 (dashed, red lines), d = 4 (dotted, blue lines), and d = 6 (dash-dot, black lines) times of the Si 
interlayer distance (3.15 A). 

The tunneling through metal|semiconductor!metal junctions can be understood in terms of the complex band 
structure of the semiconductor barrieri^ The tunneling current is carried out by evanescent states near the 
metaljsemiconductor interfaces on the semiconductor side. The energy dispersions of the evanescent and propa¬ 
gating states are called the complex band structure. The junction is assumed to be periodic in the x-y plane, so the 
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in-plane components of the fc-vector k|| = {kx,ky) are always real, while the component along the z-direction kz is 
allowed to be complex. The Fermi energy of the junction lies in the band gap region of the semiconductor barrier, so 
the imaginary part of kz is always nonzero. Although there are infinitely many kz for each k||, we only considered 
the one with the smallest imaginary part (denoted as k hereafter), since the corresponding evanescent state has the 
slowest decay rate. 

We calculated k for each ky, where all of the ky’s form a uniform mesh in the two-dimensional first Brillouin zone. 
Following Ref. [13, the density of k is defined as n(/t') = (5[K(k||) — k']. In practice the (5-function is replaced by a 

Gaussian function. The calculated density of k for AA-Si and ABC-Si are plotted in Fig. [51(a) and (b), respectively. 
The density of k in the high-K region (k > 0.15/A) is much larger than that in the Iow-k region {k < 0.15/A). 

Not all of the k’s are important for the tunneling. One k is important if the number of k is large, and/or if 
the corresponding transmission probability is large; the relative importance of each k is defined as I{k) = 

n(K)e“'''^/[n(Kniin)e“'^‘”‘“'^], where d is the thickness of the semiconductor barrierThe relative importance of the 
smallest k is scaled to be 1. For very large d the relative importance of all k except Kmin is zero, since /(k) oc 
g- 2 (K-K„in)(i^ and the tunneling is only contributed by the evanescent state corresponding to Kmin, but for thin 
barriers, k’s larger than Kmin can also be important. 

The calculated relative importance curves for several different thicknesses of the silicon barrier with AA-Si and ABC- 
Si stacking orders are shown in Fig. jSKc) and (d) respectively, for the thicknesses d of Si barriers are corresponding to 
2, 4, and 6 atomic layers. The thickness d is estimated to be smaller by 2 atomic layers than the actual thickness of the 
Si barrier, since the interface-induced changes in the potential are confined in the hrst Si atomic layer at the interface. 
Thus Figs. |SJc,d) are corresponding to junctions with 4, 6, and 8 Si atomic layers. The range of important k’s, i.e., 
the k’s with a significant /(k), is larger for AA-Si than ABC-Si, which explains the faster decay of the transmission 
of AA-Si based junctions. 


B. Resistance of junctions 
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FIG. 6. (Color online) Resistance-area products of (a) Ag|AA-Si|Ag and (b) Ag|ABC-Si|Ag junctions as a function of the 
number of Si layers calculated using the SCB equation (black rectangles) converge to those calculated using the Landauer 
formula (red dots). 


Using the Landauer formula, the resistance-area products (RA) of these junctions is 


RA = 


1 2'!Th 


A, 


( 8 ) 


where T is the averaged transmission defined in Eq. ([7]), 27r?i/2e^ = 12.9kll is the quantum resistance, and A is 
the cross section of the junction unit cell in the x-y plane. The resistance calculated using the Landauer formula 
Eq. ([8]) corresponds to the resistance measured using the “two-probe” configuration. The two-probe resistance can 
be interpreted as the corresponding four-probe resistance in series with a contact resistance. The contact resistance 
is the result of the mismatch in the number of channels in the leads and in the reservoirs: there are a finite number 
of channels in the leads, but infinite in the reservoirsi^ The two-dimensional multilayer-silicene junctions in question 
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actually have infinite channels in their leads. As a result, the Landauer formula Eq. ([8]) becomes invalid in the 
large-transmission region. 

To resolve this problem, we then used the SCB theory to calculate the RA of junctions, which correspondis to the 
resistance measured in the “four-probe” configuration. The values of RA calculated using the Landauer formula and 
the SCB theory are presented in Fig. [6l The SCB theory gives almost the same RA as the Landauer formula for 
junctions with more than five layers of silicene, where the averaged transmission per channel is about 0.05. In this 
small transmission limit, the contact resistance is negligible with respect to the junction resistance; thus the two-probe 
resistance is nearly equal to the four-probe resistance. This is the first time that the SCB theory is applied down 
to the small transmission region. We note that in this small transmission region, the Boltzmann equation becomes 
numerically difficult to solve; as a result, it is more convenient to use the Landauer formula instead. 

The ABC-Si junction with two Si atomic layers has a RA 50% smaller than that of the corresponding AA-Si 
junction, indicating that the calculated resistance using the SCB theory is not entirely determined by the averaged 
transmission. The existing various multi-channel extensions^i^S of the four-probe Landauer formula have the common 
feature that the group velocities of Bloch states in the leads play a role. The factors determining the SCB resistance are 
however rather difficult to analyze, due to the self-consistent nature of the Boltzmann equation. The most important 
observation from Fig. [S] is that the upper limit of the transmission per channel for the applicability of the Landauer 
formula is equal to 0.05. In the next section, we studied /i-BN-based junctions to explore whether this value (0.05) is 
universal. 


IV. MULTILAYER HEXAGONAL BORON NITRIDE JUNCTIONS 

In this section, we turn to junctions based on the wide-gap insulator hexagonal boron nitride (/i-BN). Hexagonal 
boron nitride is generally recognized as a good insulator because of its large energy band gap. However, first-principles 
calculations in the literature^ML showed that the transmission of monolayer /i-BN-based junctions is on the order of 
unity. According to the preceding results on multilayer silicene junctions, the Landauer formula is no longer applicable, 
but instead SCB theory is needed to calculate the resistivity of monolayer /i-BN junctions. 
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FIG. 7. (Color online) Resistance-area products of Ni(lll)|/i-BN|Ni junctions with different numbers of /i-BN layers calculated 
using the Landauer formula (black rectangles) and using the SCB theory (red dots). The averaged transmissions are denoted 
by blue circles. 

The electrode for /i-BN junctions is chosen to be /cc-Ni because of the small lattice mismatch, 2.504 A for /i-BN 
versus 2.492 A for the (lll)-surface of Ni. The in-plane lattice constant of /i-BN is squeezed slightly to match that 
of Ni(III), in order to simulate the commensurate /i-BN/Ni(Ill) interface. The distance between /i-BN layers is set 
to the interlayer distance in bulk /i-BN, 3.33 A. At the /i-BN/Ni(lIl) interface, N atoms sit on the top of surface Ni 
atoms, and the distance between /i-BN and Ni(lll) surfaces is 2.1 A4^ 

The transmissions of mono-, bi-, and tri-layer /i-BN based junctions are calculated using the same method for 
Ag I silicene I Ag junctions, except that spin-polarized calculations were carried out; the magnetic moments in the two 
Ni(lll) leads are set to be parallel to each other. The calculated transmission decays exponentially as a function of 
number of /i-BN layers (Fig. [T]), in accord with previous calculationsNo spin flip is considered in calculating the 
resistance, and the two spin channels are considered to be independent. The RA of monolayer /i-BN-based junctions 





calculated using the Landauer formula (Eq. ([8])) is about two times larger than that using the SCB theory, as shown 
in Fig. [T] The difference in RA between Landauer formula and the SCB theory becomes very small for junctions with 
more than one /i-BN layer. The junction with bilayer /i-BN, which is the threshold for the Landauer formula to be 
applicable, has a transmission of about 0.06. We note that this value is very close to the value 0.05 of the critical 
transmission for multilayer-silicene junctions. 


V. SUMMARY 

In this work we extended our previous worb22 to multilayer silicene junctions with up to eight layers, and also 
compared with /i-BN junctions. Two stacking orders AA and ABC of multilayer-silicene were considered. The 
calculated transmission decays as a function of the barrier thickness for junctions with more than four silicene layers. 
The electrical resistances were calculated using the Landauer formula and the SCB theory. We observed that the 
SCB resistance is not entirely determined by the averaged transmission. 

Most importantly, we learned from calculations on multilayer-silicene junctions that the upper limit of the transmis¬ 
sion per channel for the applicability of the Landauer formula is 0.05, above which the Landauer formula significantly 
overestimates the SCB resistance. Additional calculations on the /i-BN junctions also give a very close value (0.06) 
for the critical transmission. 

The SCB theory to calculate the four-probe resistance of junctions is revisited. We focused on the CPP configuration, 
in particular we present the numerical method to solve the Boltzmann equation in details. We believe that this work 
and our former work (Ref. [23) boost the applications of the SCB theory. 
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Appendix A: Numerical solution for the Boltzmann equation 


In the appendix the numerical method for solving the Boltzmann equation is presented. The Boltzmann equation 
for the CPP geometry is. 


d 1 




= -e?^z(k||)fz. 


(Al) 


The solution of Eq. (lAll) can be written as the sum + wU The first term satisfies Eq. (lAll) for ii{z) = 0 

d 1 


for which we can write an analytical solution for , 




M'^ (z,k||) = —eu^(k||)£’ 2 T-I-A'^ (k||) exp 


vi{k\\)i 


(A2) 


(A3) 


the are parameters to be determined by the boundary conditions. Then, obeys 

. 9w^'(2.k||) , w^(z,k||)-/r(z) 

az +-;-= "■ 


(A4) 


Here, approximating the z dependence of fj.{z) to second order. 
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we find the solution for is, 


with 


{z + Az, k ||) = (z, k||) exp 


Az 


ui(k|i)i 


- &oMo - bi^iAz - - 62^*2Az^, 


60 = exp 


Az 


ui(k||)T 


- 1 , 


6 , = _ 2 d|Ai„ _ 1 . 


(A6) 


(A7) 

(A8) 

(A9) 


The z-axis is uniformly discretized in practical calculations. The value of on a grid point can be derived from the 
value on one of its neighbors, since Eq. (IMl) is a recursive type equation. 

The boundary conditions for the distribution functions in leads are determined by the transmission and reflection 
coefficients of the junctions between leads. Here we only considered the systems consisting of one single junction 
sandwiched by two leads; the extension to multiply leads and junctions is straightforward. Suppose the junction 
between the leads is located at z = Zq. The leads above (z > zq) and below (z < zq) the junction are denoted as the 
“right” {R) and the “left” (L) leads, respectively. The boundary conditions for the distribution function are, 


Nr 

E 

j' 

Nl 


Nl 


P (^o") + E (^0 ) 

(AlO) 

J 

Nr 

2 h^/ (Zq-) + g |2 (Zo+) 

(All) 


The k|| dependence of the distribution function h and the transmission and reflection coefficients {t and r) are omitted 
in Eqs. (jAlOllATT]) . The total number of channels in the left and the right leads are 2Nl and 2Nii respectively, half 
of them propagating towards (<) the junction and the other half against (>) the junction. 

There is indeed a freedom to choose the boundary conditions for and , and in our implementation both and 
satisfy the boundary conditions Eqs. (lAlOllAlT]) . Note that in Eq. (lA.SI) is decoupled for different k||’s, so do the 
boundary conditions in Eqs. (lAlOirATTI) . At each k||, there are 2Nl and 27Vfl unknowns A-^’s in the left and right leads 
respectively, so the total number of unknowns is 2Nl + 2Nfi. The boundary conditions in Eqs. (lAlOilATTI) provide 
{Nl + Nr) equations, which is less than the number of unknowns by {Nr + Nr). In fact, we considered that the leads 
have finite lengths, ie., the left lead extends down to z = zl and the right lead up to z = zr {zr > zq > zl). The 
boundaries at z = zl and zr provide boundary conditions similar to Eqs. (lAlOilATTI) but without the transmission 
part; thus they provide another Nr + Nr equations for determining the unknowns . The reflection coefficients of 
the boundaries at z = zl and zr are arbitrary and they do not change the resistance of the junction if the length of 
both leads are long enough. 

The in Eq. (IA4I1 is no longer decoupled for different ky, because the chemical potential ^ contains a summation 
over all the ky’s. As a result, Eq. (IA4|) for is a self-consistent field equation, which can be solved using an iterative 
algorithm: given an initial guess for the chemical potential /rin(z), can be solved using the same method to solve 

. After solving for w^{z), a new chemical potential fJ.out{z) can be built and used for the next iteration step. 
The iteration stops when /iin(-z) agrees with ^Xout{z) within some predefined numerical accuracy. Efficient mixing 
algorithms such as Ref. [4^ can be used to accelerate the convergence. The typical number of iteration steps requried 
for current-density conservation is several hundreds. 


Appendix B: The expression for the local chemical potential 


We derive in this Appendix the expression for the local chemical potential /i(z). The number of electrons at position 
z with energy E is ky ky, E). The electrons at z are in local equilibrium with a local chemical potential ft(r), 

i.e., the occupation number of each mode is equal to -|- 1 ), where P = l/ksT, so we have 


E g/3[£:-At(z)] _|_ X ’ 

J,k|| j,k|| 


(Bl) 





















or, averaged over j and ky, 
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(/'(*. k|i. E)>,,k, = +1 = 

The equilibrium distribution function /o denotes electrons that are in global equilibrium with a global chemical 
potential ^oi 


ME) = TiE - ^io)- 


(B3) 


From Eqs. 


in Eq. o, we have, 


- {h^{z,kM, W„ - f,{z)] - T[E - M- 


dE 

In the linear-response regime, the value of /j,(z) remains close to p,o, and 

dE{x) 


E\E - /r(z)] - T\E - ^o] = -[d-(.z) - Mo] 


dx 


(B4) 


(B5) 


x—E—fio 


At low temperatures, we have —dfo{E)/dE = — dE{x)/dx\^^^_^^ = S{E — mo)- As a result, the expression for the 
local chemical potential is 


m ( z ) = +do- 


(B6) 


Suppose that W(z, ky) is the solution of the Boltzmann equation with d{z)- If the local chemical potential d{z) is 
shifted by a constant C, the corresponding solution of the Boltzmann equation becomes h^{z,'k^<^) -I- C. Both the 
voltage drop and the current density are independent of the constant C according to their definitions. So, we drop 
the term mo in the expression of fJ-iz) and obtain finally 

m(2 ) = (h^(2,k||))^.^^^ . (B7) 


* Corr. author: Hai-Ping Cheng, hping@ufl.edu 
^ K. Takeda and K. Shiraishi, Phys. Rev. B 50, 14916 (1994) 

^ E. Durgun, S. Tongay, and S. Ciraci, Phys. Rev. B 72 , 075420 (2005) 

® G. G. Guzman-Verri and L. G. Lew Yan Voon, Phys. Rev. B 76 , 075131 (2007) 

S. Cahangirov, M. Topsakal, E. Aktiirk, H. §ahin, and S. Giraci, Phys. Rev. Lett. 102, 236804 (2009) 

® P. Vogt, P. De Padova, C. Quaresima, J. Avila, E. Frantzeskakis, M. G. Asensio, A. Resta, B. Ealet, and G. Le Lay, 

Phys. Rev. Lett. 108, 155501 (2012) 

® L. Ghen, G.-G. Liu, B. Feng, X. He, P. Gheng, Z. Ding, S. Meng, Y. Yao, and K. Wu, Phys. Rev. Lett. 109 , 056804 (2012) 

^ A. Fleurence, R. Friedlein, T. Ozaki, H. Kawai, Y. Wang, and Y. Yamada-Takamura, Phys. Rev. Lett. 108 , 245501 (2012) 

® L. Meng, Y. Wang, L. Zhang, S. Du, R. Wu, L. Li, Y. Zhang, G. Li, H. Zhou, W. A. Hofer, and H.-J. Gao, Nano Letters 
13 , 685 (2013). 

® C.-L. Lin, R. Arafune, K. Kawahara, M. Kanno, N. Tsukahara, E. Minamitani, Y. Kim, M. Kawai, and N. Takagi, 
Phys. Rev. Lett. 110 , 076801 (2013) 

Z.-X. Guo, S. Furuya, J.-i. Iwata, and A. Oshiyama, Phys. Rev. B 87, 235435 (2013) 

Z.-X. Guo, S. Furuya, J.-i. Iwata, and A. Oshiyama, Journal of the Physical Society of Japan 82 , 063714 (2013) 

Y.-P. Wang and H.-P. Cheng, Phys. Rev. B 87 , 245430 (2013) 

P. Gori, O. Pulci, F. Ronci, S. Colonna, and F. Bechstedt, J. Appl. Phys. 114 , 113710 (2013). 

D. Tsoutsou, E. Xenogiannopoulou, E. Golias, P. Tsipas, and A. Dimoulas, Appl. Phys. Lett. 103 , 231604 (2013). 

P. De Padova, P. Vogt, A. Resta, J. Avila, 1. Razado-Golambo, G. Quaresima, C. Ottaviani, B. Olivieri, T. Bruhn, T. Hira- 
hara, T. Shirai, S. Hasegawa, M. Carmen Asensio, and G. Le Lay, Appl. Phys. Lett. 102, 163106 (2013). 

P. Vogt, P. Capiod, M. Berthe, A. Resta, P. De Padova, T. Bruhn, G. Le Lay, and B. Grandidier, Appl. Phys. Lett. 104 , 
021602 (2014). 

P. D. Padova, J. Avila, A. Resta, 1. Razado-Colambo, C. Quaresima, C. Ottaviani, B. Olivieri, T. Bruhn, P. Vogt, M. C. 
Asensio, and G. L. Lay, J. Phys.: Condens. Matter 25, 382202 (2013). 

E. Salomon, R. E. Ajjouri, G. L. Lay, and T. Angot, J. Phys.: Condens. Matter 26, 185003 (2014). 

C. Kamal, A. Chakrabarti, A. Banerjee, and S. K. Deb, J. Phys.: Condens. Matter 25, 085508 (2013). 








11 


Z.-X. Guo and A. Oshiyama, Phys. Rev. B 89, 155418 (2014) 

T. Sliirai, T. Shirasawa, T. Hiraliara, N. Fukui, T. Takaliashi, and S. Hasegawa, Phys. Rev. B 89, 241403 (2014) 

Y.-P. Wang, J. N. Fry, and H.-P. Cheng, Phys. Rev. B 88, 125428 (2013) 

P. E. Blochl, Phys. Rev. B 50, 17953 (1994) 

G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999) 

G. Kresse and J. Furthmiiller, Comput. Mater. Sci. 6, 15 (1996). 

G. Kresse and J. Furtlirniiller, Phys. Rev. B 54, 11169 (1996) 

J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996) 

S. Datta, Electronic Transport in Mesoscopic Systems (Cambridge University Press, Cambridge, 1995). 

J. Taylor, H. Guo, and J. Wang, Phys. Rev. B 63, 245407 (2001) 

Y. Xue, S. Datta, and M. A. Ratner, Chem. Phys. 281, 151 (2002). 

M. Brandbyge, J.-L. Mozos, P. Ordejon, J. Taylor, and K. Stokbro, Phys. Rev. B 65, 165401 (2002) 

J. M. Soler, E. Artacho, J. D. Gale, A. Garca, J. Junquera, P. Ordejon, and D. Snchez-Portal, 
Journal of Physics: Condensed Matter 14, 2745 (2002) 

N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991) 

W. H. Butler, X.-G. Zhang, and J. M. MacLaren, J. Superconduct. 13, 221 (2000). 

W. H. Butler, X.-G. Zhang, and J. M. MacLaren, J. Appl. Phys. 87, 5173 (2000). 

K. Umemoto, R. M. Wentzcovitch, S. Saito, and T. Miyake, Phys. Rev. Lett. 104, 125504 (2010) 

P. Mavropoulos, N. Papanikolaou, and P. H. Dederichs, Phys. Rev. Lett. 85, 1088 (2000) 

M. Biittiker, Y. Imry, R. Landauer, and S. Pinhas, Phys. Rev. B 31, 6207 (1985) 

M. D. Ventra, Eloectrical Transport in Nanoscale Systems (Cambridge University Press, 2008). 

O. V. Yazyev and A. Pasquarello, Phys. Rev. B 80, 035408 (2009) 

V. M. Karpan, P. A. Khomyakov, G. Giovannetti, A. A. Starikov, and P. J. Kelly, Phys. Rev. B 84, 153406 (2011) 

E. Rokuta, Y. Hasegawa, K. Suzuki, Y. Gamou, C. Oshima, and A. Nagashima, Phys. Rev. Lett. 79, 4609 (1997) 

D. D. Johnson, Phys. Rev. B 38, 12807 (1988) 




